LSTM AutoEncoder¶
気象情報をセンサデータと見立てて異常検知してみる。入力は気象庁の過去の天気情報で入手。1980年から2020年(06月まで)の毎日の天気情報。
In [1]:
import pandas as pd
import numpy as np
import os
from tqdm import tqdm
import tensorflow as tf
import pandas as pd
pd.options.mode.chained_assignment = None
import seaborn as sns
import matplotlib.pyplot as plt
from matplotlib.pylab import rcParams
import plotly.express as px
import plotly.graph_objects as go
%matplotlib inline
sns.set(style='whitegrid',palette ='muted')
rcParams['figure.figsize']=14, 8
np.random.seed(1)
# tf.random.set_seed(1)
print("TensorFlow Version: ",tf.__version__)
データの準備¶
予め超かんたんなクレンジングを実行済み。input配下にCSVが入れてある状況。
サンプルとして一番初めの月だけ見てみる。その後全期間のデータをマージ、風向きデータのOneHotが面倒なので可視化後に除外。数値データだけ残して正規化する。
※最高気温、最低気温、平均気温など相関の強いカラムが多いので本来はカラムを絞ったり次元削減するべきだが、とりあえずのモデルなのでそこまではこだわらない。
In [2]:
INPUT_DIR = 'input/'
In [3]:
csv_files = os.listdir(INPUT_DIR)
csv_files.sort()
In [4]:
# とりあえず初めの月だけ表示してみる
input_path = os.path.join(INPUT_DIR, csv_files[0])
df = pd.read_csv(input_path)
df
Out[4]:
In [5]:
# 全データを結合
df = pd.DataFrame(columns=df.columns)
for csv_file in tqdm(csv_files):
input_path = os.path.join(INPUT_DIR, csv_file)
tmp_df = pd.read_csv(input_path, parse_dates=['date'])
df = pd.concat([df, tmp_df])
df = df.fillna('').reset_index(drop=True)
In [6]:
# 風向きはいらない
direction_cols = ['wind_velocity_max_direction', 'maximum_instantaneous_wind_speed_direction', 'most_common_wind_direction']
train_cols = [col for col in df.columns.tolist() if not (col in direction_cols)]
df = df[train_cols]
In [7]:
train_size = int(len(df)*0.85)
test_size = len(df)-train_size
train,test = df.iloc[0:train_size],df.iloc[train_size:len(df)]
print(train.shape,test.shape)
feature_cols = train.drop('date', axis=1).columns.tolist()
可視化¶
超大雑把に可視化しておく。残念ながら積雪や降雪の記録は取れていないらしい。三浦地方は大雪は少ないから異常がわかりやすいと思ったのに。。。
pandas-profilingがものすごくリッチに成長していてビックリ。先述の通り相関が高いものを警告したりしてくれてる。
In [8]:
for col in feature_cols:
fig = go.Figure()
fig.add_trace(go.Scatter(x=df['date'],y=df[col],mode='lines',name=col))
fig.update_layout(showlegend=True)
fig.show()
In [9]:
import pandas_profiling as pdp
pdp.ProfileReport(df)
Out[9]:
正規化¶
標準化、正規化は悩みどころだけど、あとで各カラムの予測値と実測値の差分を足し合わせる関係で正規化にしておく
In [10]:
from sklearn.preprocessing import StandardScaler, MinMaxScaler
X_train = train[feature_cols]
X_test = test[feature_cols]
scaler = MinMaxScaler()
scaler = scaler.fit(X_train)
X_train = scaler.fit_transform(X_train)
X_test = scaler.fit_transform(X_test)
X_train = X_train.reshape(X_train.shape[0], 1, X_train.shape[1])
X_test = X_test.reshape(X_test.shape[0], 1, X_test.shape[1])
モデル生成¶
あまり深くする必要もないのだが、せっかくなので多少は積んでおく。推論が遅くて使えなかったら考える。
In [11]:
from tensorflow.keras.models import Sequential, Model
from tensorflow.keras.layers import Input,LSTM,Dense,Dropout,RepeatVector,TimeDistributed
from tensorflow.keras import regularizers
def build_lstm_ae(X_train):
inputs = Input(shape=(X_train.shape[1], X_train.shape[2]))
x = LSTM(256, activation='relu', return_sequences=True,
kernel_regularizer=regularizers.l2(0.00))(inputs)
x = Dropout(0.1)(x)
x = LSTM(128, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(64, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(32, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(16, activation='relu', return_sequences=False)(x)
x = Dropout(0.1)(x)
x = RepeatVector(X_train.shape[1])(x)
x = LSTM(16, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(32, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(64, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(128, activation='relu', return_sequences=True)(x)
x = Dropout(0.1)(x)
x = LSTM(256, activation='relu', return_sequences=True)(x)
output = TimeDistributed(Dense(X_train.shape[2]))(x)
model = Model(inputs=inputs, outputs=output)
return model
model = build_lstm_ae(X_train)
model.compile(loss='mae',optimizer='adam',metrics=['accuracy'])
model.summary()
In [12]:
# AuroEncoderなのでラベルは不要。特徴量とラベルがおんなじって感じ。仮にノイズ除去モデルを作りたいなら第一引数をノイズを付与したX_trainにする
es = tf.keras.callbacks.EarlyStopping(monitor='val_loss',patience=5,mode='min')
history = model.fit(
X_train, X_train, epochs=100,
batch_size=256, validation_split=0.1, callbacks=[es], shuffle=False
)
In [13]:
plt.plot(history.history['loss'],label='Training Loss')
plt.plot(history.history['val_loss'],label='Validation Loss')
plt.legend()
Out[13]:
In [14]:
model.evaluate(X_test,X_test)
Out[14]:
異常検知¶
AutoEncoderなので予測値は実際の気象情報を変換して復元した値になる。復元した値と実際の値の差分を取って、差分が大きいところを異常とする。
気象情報なので多少は異常があると仮定しているが、仮に本物のセンサデータで異常発生時のデータが無い場合は差分の平均との差分が+3σより大きかったら異常とか、観測データの差分の最大値+αを超えたら異常って扱いにするのかな。
差分はMAEとMSEで見ているが、正規化しているから差分は1未満になるのでMAEよりMSEのほうが顕著に出ているように見える。
In [15]:
# plot the loss distribution of the training set
X_pred = model.predict(X_train)
X_pred = X_pred.reshape(X_pred.shape[0], X_pred.shape[2])
X_pred = pd.DataFrame(X_pred, columns=feature_cols)
X_pred.index = train.index
scored = pd.DataFrame(index=train.index)
Xtrain = X_train.reshape(X_train.shape[0], X_train.shape[2])
scored['Loss_mae'] = np.mean(np.abs(X_pred-Xtrain), axis = 1)
scored['Loss_mse'] = np.mean(np.power(X_pred-Xtrain, 2), axis=1)
plt.figure(figsize=(16,9), dpi=80)
plt.title('Loss Distribution', fontsize=16)
sns.distplot(scored['Loss_mae'], bins=50, kde=True, color='red', label='mae')
sns.distplot(scored['Loss_mse'], bins=50, kde=True, color='blue', label='mse')
plt.xlim([0.0,.3])
Out[15]:
In [16]:
X_pred = model.predict(X_test)
X_pred = X_pred.reshape(X_pred.shape[0], X_pred.shape[2])
X_pred = pd.DataFrame(X_pred, columns=feature_cols)
X_pred.index = test.index
scored = pd.DataFrame(index=test.index)
Xtest = X_test.reshape(X_test.shape[0], X_test.shape[2])
scored['mae'] = np.mean(np.abs(X_pred-Xtest), axis = 1)
scored['mse'] = np.mean(np.power(X_pred-Xtest, 2), axis = 1)
plt.figure(figsize=(16,9), dpi=80)
plt.title('Test Error Distribution', fontsize=16)
sns.distplot(scored['mae'], bins=50, kde=True, color='red', label='mae')
sns.distplot(scored['mse'], bins=50, kde=True, color='blue', label='mse')
plt.xlim([0.0,.3])
Out[16]:
In [17]:
X_test_pred = model.predict(X_test)
mse = np.mean(np.power(X_test - X_test_pred, 2), axis=1)
mse = np.sum(mse, axis=1)
mae = np.mean(np.abs(X_test - X_test_pred), axis = 1)
mae = np.sum(mae, axis=1)
test_df = test.copy()
test_df['mse'] = mse
test_df['mae'] = mae
In [18]:
# 全体のthreshold*100 % は正常値として扱い、それ以外を異常とする。そのための閾値を計算する。
threshold = 0.998
mse_threshold = np.quantile(mse, threshold)
mae_threshold = np.quantile(mae, threshold)
print('MSE {0} threshold:{1}'.format(threshold, mse_threshold))
print('MAE {0} threshold:{1}'.format(threshold, mae_threshold))
In [19]:
display(test_df[test_df['mse'] > mse_threshold])
display(test_df[test_df['mae'] > mae_threshold])
In [20]:
test_df['mse_threshold'] = mse_threshold
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse'],mode='lines',name='mse'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse_threshold'],mode='lines',name='mse_threshold'))
fig.update_layout(showlegend=True)
fig.show()
test_df['mae_threshold'] = mae_threshold
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae'],mode='lines',name='mae'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae_threshold'],mode='lines',name='mae_threshold'))
fig.update_layout(showlegend=True)
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse'],mode='lines',name='mse'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mse_threshold'],mode='lines',name='mse_threshold'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae'],mode='lines',name='mae'))
fig.add_trace(go.Scatter(x=test_df['date'],y=test_df['mae_threshold'],mode='lines',name='mae_threshold'))
fig.update_layout(showlegend=True)
fig.show()